Preprocessing report


# Packages ----------------------------------------------------------------

# using a reproducible environment
renv::restore()

# pacman allows to check/install/load packages with a single call
# if (!require("pacman")) install.packages("pacman") # already in renv.lock
library("pacman")

# packages to load (and install if needed) -------------------------------
pacman::p_load(
  here,      # easy file paths
  see,       # theme_modern and okabeito palette
  report,    # reporting various info 
  labelled,  # labelled data
  # ---- packages specific to this project ----
  scales,
  lme4,
  simr,
  readxl,
  openxlsx,
  statmod,
  qqplotr,
  emmeans,
  ggbeeswarm,
  ggpubr,
  patchwork,
  quarto,
  easystats,
  # ---
  tidyverse  # modern R ecosystem
)


# Custom functions shared across scripts ----------------------------------
source(here("scripts/_functions.R"))


# Global cosmetic theme ---------------------------------------------------

theme_set(theme_modern(base_size = 14)) # from see in easystats

# setting my favourite palettes as ggplot2 defaults
options( 
  ggplot2.discrete.colour   = scale_colour_okabeito,
  ggplot2.discrete.fill     = scale_fill_okabeito,
  ggplot2.continuous.colour = scale_colour_viridis_c,
  ggplot2.continuous.fill   = scale_fill_viridis_c
)


# Fixing a seed for reproducibility ---------------------------------------
set.seed(14051998)


# Adding all packages' citations to a .bib --------------------------------
knitr::write_bib(c(.packages()), file = here("bibliography/packages.bib"))

Down below is the code for the outlier selection procedure described in the manuscript.

Preparing the data for the outlier detection procedure

# ═══ Accuracy outliers analysis ═══════════════════════════════════════════════

df_0e <- 
  df_explicit |> 
  mutate(across(c(aphantasia, color, congruence), as.factor))

df_0i <- df_implicit |> 
  mutate(across(c(aphantasia, color, congruence), as.factor))

Listing error rates per subject

# ─── Listing error rates per subject ───
list(
  df_0e = df_0e,
  df_0i = df_0i
  ) |>
  imap(
    ~.x |>
      select(subjectid, starts_with("correct"), aphantasia) |>
      group_by(subjectid) |>
      count(pick(2)) |>
      filter(pick(1) == 1) |>
      ungroup() |>
      mutate(
        n_tot = max(n),
        prop = (n_tot - n) / n_tot * 100  # analysing the error rate per subject
      ) %>%
      arrange(desc(prop)) |>
      select(1, 5)
    )
# 5 accuracy outliers in the explicit task, 4 in the implicit one

Removing error trials and irrelevant variables

# ─── Removing incorrect trials and task-wise accuracy outliers ───
df_1e_out <-
  df_0e |> 
  # filtering out...
  filter(
    # participants identified with with high error rates
    !(subjectid %in% c( 
      "subject_7",
      "subject_94", 
      "subject_25", 
      "subject_4",
      "subject_97"))
    ) |>  
  # removing irrelevant variables
  select(-c(sex, vviq80, orientation, response))

df_1i_out <-
  df_0i |> 
  filter(
    !(subjectid %in% c(
      "subject_21",
      "subject_56",
      "subject_9"
    ))
  ) |>  
  # removing irrelevant variables
  select(-c(sex, vviq80, orientation, response))

df_1e <-
  df_1e_out |> 
  # filtering out...
  filter(correct_explicit == 1) |> 
  select(!correct_explicit)

df_1i <-
  df_1i_out |> 
  # filtering out...
  filter(correct_implicit == 1) |> 
  select(!correct_implicit)

Removing extreme RTs

# ═══ Broad RT outlier trials removal ════════════════════════════════════

df_2e <- 
  df_1e |>  
  # filtering out extreme RTs
  filter(rt > .3 & rt < 3)

df_2i <- 
  df_1i |>  
  # filtering out extreme RTs
  filter(rt > .3 & rt < 3)

Creating a table with RT means

# ═══ RT means table - Both tasks ══════════════════════════════════════════════

df_rt <- 
  list(df_2e = df_2e, df_2i = df_2i) |> 
  imap(
    ~.x |> 
      select(rt) |> 
      report() |>
      as.data.frame() |>  
      select(1:10)
  )

Detecting outliers in the RT means

# ═══ Aberrant RT means outliers removal ═══════════════════════════════════════

# ─── Finding outliers ───
threshold_e <- df_rt$df_2e$Median + 3*df_rt$df_2e$MAD
threshold_i <- df_rt$df_2i$Median + 3*df_rt$df_2i$MAD

df_2e |>
  group_by(subjectid) |>
  summarise(mean_rt = mean(rt)) |>
  ungroup() |>
  filter(mean_rt > threshold_e) |>
  arrange(desc(mean_rt))
# 7 outliers in the explicit task

df_2i |>
  group_by(subjectid) |>
  summarise(mean_rt = mean(rt)) |>
  ungroup() |>
  filter(mean_rt > threshold_i) |>
  arrange(desc(mean_rt))

Removing outliers in the RT means

# ─── Removing specific outliers ───

df_3e <- 
  df_2e |> 
  filter(
    !(subjectid %in% c(
      "subject_49",
      "subject_59",
      "subject_107",
      "subject_100",
      "subject_73",
      "subject_106",
      "subject_119"
    ))
  )

df_3i <- 
  df_2i |> 
  filter(
    !(subjectid %in% c(
      "subject_49",
      "subject_107",
      "subject_30",
      "subject_120",
      "subject_127"
    ))
  )

Preparing the data for the correlation analysis

# ─── Preparing tables for correlation ───
congruence_effects_emm <-
  list(
      df_3e = df_3e,
      df_3i = df_3i
      ) |> 
  imap(
    ~.x |> 
      group_by(subjectid, congruence, color) |> 
      reframe(mean_rt = mean(rt)) |> 
      group_by(subjectid, congruence) |> 
      reframe(mean = mean(mean_rt)) |> 
      pivot_wider(
        names_from = congruence,
        values_from = mean
      ) |> 
      mutate(congruence_effect = uncongruent - congruent, .keep = "unused") |> 
      left_join(df_questionnaires[,c(1,5:8)], by = "subjectid") |> 
      ungroup()
    )

Creating extreme groups

df_3ex <-
  df_3e |> 
  left_join(df_questionnaires[c(1,5)], by = "subjectid") |> 
  filter(vviq80 == 16 | vviq80 >= 46)

df_3ix <-
  df_3i |>
  left_join(df_questionnaires[c(1,5)], by = "subjectid") |> 
  filter(vviq80 == 16 | vviq80 >= 42)

df_3ix |> 
  group_by(subjectid, aphantasia) |> 
  count() |> 
  group_by(aphantasia) |> 
  count()

Session information for reproducibility

═════════════════════════════════════════════════════════════════════════

Analyses were conducted using the R Statistical language (version 4.4.1; R Core
Team, 2024) on Windows 11 x64 (build 22631)

Packages used:

  - quarto (version 1.4.4; Allaire J, Dervieux C, 2024)
  - qqplotr (version 0.0.6; Almeida A et al., 2018)
  - lme4 (version 1.1.35.5; Bates D et al., 2015)
  - Matrix (version 1.7.0; Bates D et al., 2024)
  - effectsize (version 0.8.9; Ben-Shachar MS et al., 2020)
  - ggbeeswarm (version 0.7.2; Clarke E et al., 2023)
  - statmod (version 1.5.0; Giner G, Smyth GK, 2016)
  - simr (version 1.0.7; Green P, MacLeod CJ, 2016)
  - lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  - ggpubr (version 0.6.0; Kassambara A, 2023)
  - labelled (version 2.13.0; Larmarange J, 2024)
  - emmeans (version 1.10.3; Lenth R, 2024)
  - parameters (version 0.22.0; Lüdecke D et al., 2020)
  - performance (version 0.12.0; Lüdecke D et al., 2021)
  - easystats (version 0.7.2; Lüdecke D et al., 2022)
  - see (version 0.8.4; Lüdecke D et al., 2021)
  - insight (version 0.20.1; Lüdecke D et al., 2019)
  - bayestestR (version 0.13.2; Makowski D et al., 2019)
  - modelbased (version 0.8.8; Makowski D et al., 2020)
  - report (version 0.5.8; Makowski D et al., 2023)
  - correlation (version 0.8.5; Makowski D et al., 2022)
  - here (version 1.0.1; Müller K, 2020)
  - tibble (version 3.2.1; Müller K, Wickham H, 2023)
  - datawizard (version 0.11.0; Patil I et al., 2022)
  - patchwork (version 1.2.0; Pedersen T, 2024)
  - R (version 4.4.1; R Core Team, 2024)
  - pacman (version 0.5.1; Rinker TW, Kurkiewicz D, 2018)
  - openxlsx (version 4.2.5.2; Schauberger P, Walker A, 2023)
  - ggplot2 (version 3.5.1; Wickham H, 2016)
  - forcats (version 1.0.0; Wickham H, 2023)
  - stringr (version 1.5.1; Wickham H, 2023)
  - tidyverse (version 2.0.0; Wickham H et al., 2019)
  - readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  - dplyr (version 1.1.4; Wickham H et al., 2023)
  - purrr (version 1.0.2; Wickham H, Henry L, 2023)
  - readr (version 2.1.5; Wickham H et al., 2024)
  - scales (version 1.3.0; Wickham H et al., 2023)
  - tidyr (version 1.3.1; Wickham H et al., 2024)

═════════════════════════════════════════════════════════════════════════